Multimode number-squeezing in collisions of two Bose-Einstein condensates 



P. Deuar, 1 T. Wasak, 2 P. Zin, 3,4 J. Chwcdehczuk, 2 and M. Trippcnbach 2 

1 Institute of Physics, Polish Academy of Sciences, Al. Lotnikow 32/46, 02-668 Warsaw, Poland 
^Institute of Theoretical Physics, Physics Department, 
University of Warsaw, Hoza 69, PL-00-681 Warsaw, Poland 
3 Narodowe Centrum Badan Jqdrowych, ul. Hoza 69, 00-681 Warszawa, Poland 
4 Univ. Paris Sud, CNRS, LPTMS, UMR 8626, Orsay 91405 France 

We investigate supersonic collisions of two Bose-Einstein condensates as a potential source of 
entangled atomic pairs by analyzing the reduction of the number difference fluctuations between 
regions of opposite momenta. We show that due to the non-monochromaticity of the mother clouds, 
a high pair-production rate is disadvantageous for number squeezing. In the presence of many 
competing modes, smaller cloud densities are preferable. Using a simple model, we explain the 
relationship between density correlations and the number squeezing, and argue that the multi- 
mode nature of the scattering must be taken into account to explain the properties of the halo of 
scattered atoms. We also analyze the impact of the Bose enhancement on the number squeezing, 
by introducing a simplified low-gain model. 
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I. INTRODUCTION AND OUTLINE 

A supersonic collision of two Bosc-Einstcin conden- 
sates is a source of strongly correlated atomic pairs, which 
may be potentially used to create spatially-separated en- 
tangled states of massive particles Such states 
could be used to extend the study of the Einstein- 
Podolsky-Rosen paradox [j| [l(|, local realism and 
Bell inequality tests into a regime where rest mass is non- 
negligible. In the context of quantum metrology, usefully 
entangled states allow one to surpass the Standard Quan- 
tum Limit - the maximum parameter estimation preci- 
sion allowed by classical physics [lU, [l3| . 

Scattering of atomic pairs can lead to reduced fluctua- 
tions of the relative population between two regions of op- 
posite momenta. This effect, called the number-squeezing 
and a closely related violation of the Cauchy-Schwartz 
inequality have been recently observed in experiments 
[HQ- Number squeezing, if accompanied by sufficiently 
high mutual coherence between the regions, is indicative 
of spin squeezing fbH 15:]. Spin squeezed states, which are 
known to be usefully entangled from a quantum metrol- 
ogy point of view 0, 0, EH , have been recently engi- 
neered in a number of experiments fl7l - [22| . In addition, 
the perfectly number-squeezed "twin-Fock" state, which 
is not spin-squeezed but is nevertheless strongly entan- 
gled, has been recently observed by Liicke et al [23| . 

In this theoretical study we consider in detail the 
number difference squeezing between atoms with op- 
posite momenta in a collision of two realistic non- 
monochromatic metastable 4 He BECs [H,Q. We simulate 
the collision using the positive-P method in the Bogoli- 
ubov approximation [241 ] . We also introduce a simple 
yet intuitive model, which relates the number squeezing 
to the second order correlations between the scattered 
atoms and demonstrate its validity in a wide range of 
parameters. The main conclusion of this work is that the 
non-monochromatic nature of the mother clouds is the 



main limiting factor to strong number-squeezing in the 
scattering halo. We argue that in the presence of many 
competing modes, smaller cloud densities are advanta- 
geous, because mode-mixing effects are destructive for 
the number-squeezing and they become more pronounced 
with higher cloud density. 

The manuscript is organized as follows. In Section |TT] 
we introduce a model, which describes the creation of 
pairs in BECs collisions. We discuss the relevant physical 
parameters and develop an alternative low-gain model, 
useful for simulations when the bosonic amplification of 
pair scattering can be neglected. In Section HTT1 we calcu- 
late the second order correlation function of the scattered 
atoms and explain how it consists of two parts - co- linear 



and back-to-back momentum correlations, like in [251 ] . In 
Section IIVI we calculate the number-squeezing parame- 
ter and using a simple Gaussian model relate it to the 
second-order coherence of the system. In Section [V] we 
present our numerical results in the regime of both high- 
and low-density of the mother BECs and compare these 
results with the model. In the course of this analysis, the 
factors affecting the number squeezing become apparent. 
We conclude in Section IVT1 Some technical details of the 
calculations and numerics presented in the main text are 
discussed in the Appendices. 



II. THEORETICAL MODEL FOR BEC 
COLLISION 



A single stationary BEC can be split into a superposi- 
tion of two counter-propagating wave-packets by means 
of Bragg scattering J26[. In a center-of-mass reference 
frame, the average velocity of each cloud is ±v rec - a re- 
coil velocity due to an absorption and subsequent emis- 
sion of a Bragg photon. If the relative velocity 2v Tec is 
approximately above the speed of sound c = ^Jgp/m, i.e. 
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when the Mach number 

Ma = 2v rcc /c > 1, 



-1.£ 



elastic collisions of particles from the two clouds lead o.5 
to scattering of atomic pairs out of the BECs. Here 
p = N '/ '(4/37ri?jp F ) is the average density of an isotropic 
condensate in the Thomas- Fermi approximation with ra- _^>- 

diuS i?TF- 

The dynamics of the scattering is governed by the "0.5 
energy- and momentum-conservation laws, v\ + v\ ~ 
2Vf CC and vi + V2 ~ 0, where index (1,2) labels the .•) 
pair members. The equalities are only approximate, due 
to the momentum spread of the initial BEC and finite 
duration of the collision. These conservation laws dic- 
tate that atoms are scattered onto a shell of radius w roc 
centered around zero, called the scattering halo. This 
phenomenon was observed in many experiments 0-Q|- 
Atoms are usually registered after a long time of free ex- 1 
pansion, when their positions approximately correspond 
to the momentum distribution just after the collision. p. 5 

Since particles scatter in pairs, there is an expectation 
that the measured population difference between two op- 
posite regions may fluctuate below the shot-noise limit. 
In the idealized case of scattering from two plane waves, 
these fluctuations are suppressed down to zero, in anal- -0.5 
ogy to the simplest model of parametric down-conversion. 
Our study takes on the task of generalizing this simple _■] 
picture and calculating the number-squeezing in conden- 
sate collisions assuming a realistic shape and time evo- 
lution of the source, and including the time-dependent 1 q g q q g 
interactions with the mean field after the scattering (2?} ■ k x / k 



A. Collision parameters 

To focus on the essential features of the process, we 
consider the simplest case of the initial condensate pre- 
pared in the ground state of a spherical trap. Depending 
on the non-linearity and the duration of the collision, 
the scattering of atoms can be either spontaneous or en- 
ter the Bose enhanced regime. Therefore, we introduce 
a dimcnsionlcss parameter 

7 = *col/*nl 

to distinguish between these two possible scenarios. Here, 
tool = ^tfA'i-cc is the duration of the scattering process, 
while t n \ = h/(^gp) is the rate of the two-body colli- 
sions. The one-third in the denominator approximately 
accounts for the fact that collisions between the atoms 
take place mostly in the center of the trap, where the 
density is high. It has been demonstrated [27H31J that 
when 7 > 1, the system enters the stimulated regime. 

We perform numerical simulations using parameters 
of metastable 4 He atoms, i.e. m = 6.65 x 10~ 27 kg, 
a s = 7.5 x 1CP 9 m and take v Tec = 10 cm/s. To address 
both spontaneous and stimulated regimes we consider the 
following cases. In the first, which we call the dilute gas 



FIG. 1. Outcome of a single realization of a collision of two 
BECs in the dense 7 = 1.02 case (In a simulation of (|5I4[I ). 
Shown are cross-sections through the density of atoms after 
the end of the collision (t = 1 .6t co ii ) , inihek x — k y momentum 
plane (top), and k z — k y plane (bottom). Clearly, a halo 
of scattered atoms is forming, nevertheless the atoms at the 
opposite regions of the halo are poorly aligned, which implies 
the absence of number squeezing in the system. Data scaled 
to units of [feo] -3 - The white areas are saturated at this colour 
scale, the bold (thin) black contours are at 0.1 (0.01) of the 
peak condensate density. 



case, the sample consists of N = 17 540 atoms in the ini- 
tial BEC, and u = 2n x 928 s" 1 . Here Ma = 13.12 and 
7 = 0.24, hence the scattering is spontaneous all along. 
In the opposite dense gas case, we take N = 74 360 and 
ui = 2n x 1911 s _1 , which gives Ma = 6.56 and 7 = 1.02, 
so the Bosonic enhancement becomes significant. In both 
cases, we simulate the scattering process until t = 1.7t co ii 
- a time at which the collision is completed. In Fig. [TJ 
cuts through the halo density are shown for the dense case 
at the end of the collision. The result was obtained us- 
ing the positive-P numerical method, which is discussed 
in detail below. Note that although the ensemble aver- 
age of the momentum distribution is symmetric around 
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zero, isolated density grains are present in a single re- 
alization, and they show an evident and massive lack of 
the symmetry. This is related to increased fluctuations of 
the population imbalance between regions with opposite 
momentum, and will be investigated further below. 



B. BEC wave-function 

We model the initial trapped BEC by a c-number wave- 
function <^o( x ) which is a solution of the stationary Gross- 
Pitaevskii equation, 



M0o(x) 



2m 



v 2 + y trap (x)+. g |0 o (x)p 



<Mx). (l) 



of two complex fields ?/;(x, t) and ?/>(x, t). The Bogoli- 
ubov equation ([3]) corresponds to a pair of stochastic Ito 
equations, 

ihd t ^(x, t) = (-^-V 2 + 2<7|0(x, t)\A ^(x, f) (5a) 
\ Ira ) 

+# 2 (x,<)^(x,t)* + v/i^(x,t)£(x,t),(5b) 
ihdMx,t) = +2g\<f>(pc,t)A $(x,t) (5c) 



-# 2 (x,i)V>(x,i)* 



/ i%(/)(x,i)5(x,t).(5d) 



Here £(x, t) and £(x, t) are independent, real stochas- 
tic Gaussian noise fields^ with zero mean and second mo- 
ments given by (£(x, i)£(x', t')) = and 



Here m is the atomic mass, /i is the chemical potential, 
g = 4irh 2 a s /m the interaction strength related to s-wave 
scattering length a s , and V tTap (x.) = ^niu> 2 x. 2 is the har- 
monic trapping potential with a frequency ui. The col- 
lision is triggered by a pair of brief Bragg pulses shined 
onto the condensate [26|. Consequently, a superposition 
of two counter-propagating, mutually coherent atomic 
clouds is created 



ik Z 



-ikoz 



4>{x,t = 0) = <Mx) ^= = 0+(x) + 0_(x), (2) 

where fco = mv ICC /h is the wave- vector associated with 
the recoil velocity. After the pulses are applied, the trap 
is switched off, and the two condensates begin to move 
apart, activating the collision process. 



C. Positive-P method 

To describe the scattered atoms we introduce a bosonic 
operator <5(x, t). In the spirit of the Bogoliubov approx- 
imation, we use the linearized equations of motion for 
the field S, assuming that the number of scattered atoms 
is a small fraction of the condensate population and the 
self-interaction of 5 can be neglected, 



ihd t 5(x,t) = 



-^V 2 + 2 ff |0(x,t)| 2 
^ 2 (x,t)<5t(x,i), 



<*(x,t) 



(3) 



The coherent mean field part </>(x, t) during the collision 
evolves according to the time-dependent Gross-Pitacvskii 
equation, 



ihdt<t>{x., t) 



2m 



<#(x,*)| J 



(x,t). (4) 



To solve the coupled equations (j3)) and ((4]), we use the 
stochastic positive-P method [24|, where instead of di- 
rectly solving Eq. ([3]) for 6 we sample the distribution 



(£(x, i)£(x', f )) - (£(x, i)£(x', t')) = «J(x - x')6(t - t'). 

(6) 

Within the stochastic positive-P method, any physical 
quantity can be obtained with the mapping 5 — > and 
5' — > ip* , and replacing quantum averages of normally- 
ordered operators by stochastic averages [32j . Equations 
([5]) recover the exact quantum dynamics of Eq. © in the 
limit of an infinite number of samples. 



D. Scattering in the absence of bosonic 
enhancement 

The main goal of this study is to investigate how the 
number squeezing between two regions of the halo is af- 
fected by various phenomena which occur during the col- 
lision. Various phenomena influence the dynamics in a 
complicated way and it seems very advantageous if we 
could, at least theoretically, "turn on/off" some of them 
to isolate their effects. For instance, the impact of the 
BEC mean-field on the scattered particles can be easily 
"controlled" by including or excluding the second term 
of the right-hand-sides in lines (f5"a)) and (|5c]) . Similarly, 
the mean-field repulsion in the evolution of the source 
BEC can be controlled by setting g = by hand in 
Eq. (jj). One can further simplify the dynamics by mod- 
cling the two counter-propagating condensates with non- 
expanding Gaussians, substituted for </>(x, t) into Eqs ([5]). 

In this subsection we propose a simple perturbative 
model, which describes the condensate collision in the 
absence of bosonic stimulation. We start by introducing 
a hierarchy of fields 



5(x,i)=£s«(x,t). 



(7) 



The lowest order term 5^ (x, t) is the solution of the free 
equation (i.e. without the additional particle creation 
term in the Bogoliubov field) 



ihd t S w (x,t) = J ff (x,t),5 (0) (x,t), 



(8) 
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where H (x,t) = -|^V 2 + 2g\<j>(x, t)\ 2 . The higher and M, the anomalous density [29l l3ll 



is 



terms of expansion ([7]) evolve according to 

ih d t 5 (3) (x, t) = H (x, t)S {j} (x, t) + , 9 2 (x, t)tfCi-Dt( rj t y 

(9) 

In this approach, the bosonic enhancement appears only 
in the higher order fields and can be excluded by restrict- 
ing the dynamics only to the two lowest ones, namely 
j = and j = 1. 

As we argue in [33} . the set of two coupled equa- 
tions for j = and j = 1 can be formally solved 
by replacing the operators with complex stochastic 
fields, i.e. S^(x,t) <5«(x,t) and 6^(x,t) 
<5"'(x, t)*. Then, the c-number equivalent of (|8|) and 
@ is solved numerically from the initial conditions con- 
sisting of 6^(x, 0) = and setting j(°)(x, 0) as a 
random complex Gaussian field with zero mean and 
the variances (5^(x, 0)*6^(x', 0)) = 5(x - x') and 
((5 (0) (x,0)<5 (0) (x',0)) = 0. It is important to note that 
- contrary to the positive-P equations ([5]) - the stochas- 
ticity is introduced only through the initial conditions, 
very much like in the truncated Wigner method, but with 
twice the variance. 

With the solutions <5(°)(x, t) and S^'(x, t) at hand, one 
can reproduce the observables. For instance the lowest 
order correlation function is given by 



(10) 



where the over-bar denotes averaging over the ensemble 
of initial conditions. 

In the following Section, we analyze some formal prop- 
erties of the Bogoliubov equation © and introduce the 
second order correlation function of scattered atoms, 
which, as will be argued in Section IIV1 determines the 
amount of number-squeezing in the system. 



III. TWO-BODY CORRELATIONS OF 
SCATTERED ATOMS 

The normalized second order correlation function of 
scattered atoms in momentum space is defined as 

((5t(k)J(k))(Jt(k')(J(k')> 

Henceforth, we omit the explicit time dependence from 
the operator 6(k, t). Note that since the Bogoliubov 
equation of motion ([3j> is linear and the initial state of 
scattered atoms is a vacuum, with the help of Wick's 
theorem we can write 

. W(k |GW(k,k T + |M(k,k T 

g (k,k)-l+ G(1)(k;k)G(1 ) (k/ik0 • ("J 

Here, is the one-body density matrix of the scattered 
atoms defined as 



M(k,k') = (<5(k)<5(k')>. 



(14) 



In order to find a natural interpretation of the compo- 
nents of (|12[) we make some further simplifications in our 
model, that are used only in this Section. First, consider 
the following simplified model [28|, US EI of the collision 
dynamics, described by the equation 



ihd t S(x,t) = — 



h 2 V< 

2m 



«(x,t) + 2^ + (x,t)0_(x,t)*t( x>t ) J 

(15) 

Here, a pair of atoms is taken from counter-propagating 
condensates 4>± - defined in Eq. ([2]) - and placed in the 
field of scattered atoms. As compared with Eq. ©, this 
model neglects the impact of the mean field of the collid- 
ing BECs on the scattering process, as well as terms pro- 
portional to 4>± that are strongly non energy-conserving 
in the halo. 

Now let us make a second simplification regarding the 
internal dynamics of the two colliding wave-packets. In 
general the two functions 4>±(x, t) evolve according to 
the Eq. (QJ, but for the sake of the present considera- 
tions we neglect the non-linear term and use the equa- 
tions of motion of free expansion. In this case, 4>± (k, t) = 

0±(k)e _i "2^ t . Such a "reduced Bogoliubov" model with 
these two approximations has been widely used previ- 
ously to investigate the dynamics of the collisi on |28l - 
HH, j an d was investigated in some detail in |27|. 

Equation (fT5|) can be solved up to the first order in 
the perturbative regime, to obtain at times long after 
the collision 

M(k,k') = ^|2 J y"d 3 k 1 d 3 k 2 0_(k 1 )0 + (k 2 ) x (16) 

xS(k 2 + k 2 -k 2 - k' 2 )6 {3) (ki + k a - k - k'). 

The anomalous density M (k, k') can be interpreted as 
the probability amplitude for having one atom with mo- 
mentum k and the second with k'. These two momenta 
come from a coherent superposition of probability am- 
plitudes describing elementary collision events (energy 
and momentum conservation laws are satisfied) between 
atoms from the BECs that come with probability ampli- 
tudes 4>± . Since the condensate functions <fi± in (p~6|) are 
localized around ±fc , the anomalous density M(k, k') is 
non- vanishing only when k ~ — k' and |k| ~ fco- In this 
sense, M describes scattered atomic pairs, correlated for 
opposite momenta. Related arguments have been pre- 
sented in p5| . 

Using similar arguments, we obtain a useful relation 
that is valid within Bogoliubov theory in the perturbative 
regime 



G (1) (k,k') = J d 3 k"M*(k,k")M(k",k'). 



(17) 



G«(k,k')H<5t(k)<5(k')), 



(13) 



As we argued above, the first anomalous density gives the 
contribution to the integral when k ~ — k" and the sec- 
ond when k" ~ — k', hence the one-body density matrix 
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is non-zero only if k ~ k' and |k| ~ fc . In conclusion, 
the scattered atoms are described either by the co-linear 
part for k ~ k', 



(2; 



9o\ 



(k,k')=g (2) (k,k'«k)«l 



G«(k, k)G«(k',k') 
(18) 

or the back-to-back part when their wave- vectors are al- 
most opposite k ~ — k', 



^ 2 b ) (k,k')=.g (2) (k,k'«-k)«l- 



|A/(k,k') 



G«(k, k)G(!)(k',k')' 
(19) 

In our numerical simulations we have seen that the 
above interpretation of M and G' 1 ) is valid to a high de- 
gree of accuracy for a wide range of parameters, even 
when the assumptions introduced above are not fully 
valid. Therefore, throughout the rest this work, we use 
the division of the second order correlation function into 
the co-linear (|T5|) and back-to-back part ([THl) . 

At this stage we are ready to introduce the number- 
squeezing parameter and show how it relates to the g^ 2 ' 
correlation function in various relevant regimes. 



IV. NUMBER SQUEEZING IN A MULTI-MODE 
SYSTEM 

In this section we define the number-squeezing param- 
eter and show how it is related to the second-order co- 
herence of the system. 

Atoms are registered (and counted) in two bins, say a 
and b, encompassing volumes V a /f, in momentum space. 
The corresponding atom-number operators read 



i a/b = Jd 3 k5\k)6(k). 



(20) 



We introduce the number-difference operator n = h a — hb 
and define the number-squeezing parameter as follows 



(21) 



where n = (h a ) + (hb) is the total number of particles 
in both bins. States that have sub-Poissonian popula- 
tion imbalance fluctuations rj 2 < 1 are called number 
squeezed. In the symmetric case (h a ) — (fib), num- 
ber squeezing is equivalent to violation of the Cauchy- 
Schwartz inequality Q, which implies the presence of 
non-classical correlations in the system. 

Using Eq. (|2T))) and the definition of g 1 - 2 * 1 from Eq. flTTl) . 
we obtain that 



G a a + Gbb — 2G a (, 



(22) 



where the G,,- stands for a two-fold integral 



G i:j = jd 3 kjd 3 k' g^ 2 \k,k')n(k)n(k'). (23) 



Here n(k) = (jt(k)<5(k)}. Note that if the two regions a 
and b are located on the opposite sides of the halo, G aa 

(2) 

and Gbb depend on the co- linear correlation function g^ , 

(2) 

while Gab is a functional of g bb . 

To handle the integrals (|2"3")l and evaluate the number- 
squeezing parameter rj 2 for typical situations, let us 
model the normalized co-linear and back-to-back corre- 
lation functions by Gaussians 



9k 



for k ~ k' and 



S^(k,k') = l 



(2) (k,k') = 1 



n 



hbb I 

i— z,t : r 



(fci+fc^) 2 



(24) 



(25) 



for k ~ — k'. This way we only need to extract the am- 
plitudes ft. c i/bb an d the widths af^ hh from the results of 
the numerical simulations, similarly to the analysis of the 
Cauchy- Schwartz violation inQ. The Gaussian form is a 
very good match to the calculated and also the experi- 
mental correlation shapes. 



Y=1.02 
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t [tcoll 
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FIG. 2. The number squeezing parameter rf in the dense 7 = 
1.02 case as a function of time for three different bin volumes 
given in kg units. The solid lines are from a numerical solution 
of the full model (|5l4p . The dashed lines are predictions of the 
Gaussian correlation model (1261). The dash-dotted horizontal 



line denotes the shot-noise level r\ =1. 

The product runs over three orthogonal directions 
where the axis z is along the collision direction, while 
r and t arc orthogonal to each other and lie in the x — y 
plane, corresponding to radial (r) and tangential (t) di- 
rections. Additionally, we model the density n(k) in (|2"3")l 
in the following way. We assume that the bin widths 
L z and L t (in the z and t directions respectively) are 
small enough for the density to be practically constant. 
This assumption is in our case well satisfied. On the 
other hand, the density quickly decays in the r direction. 
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To account for this drop relatively simply, we model the 
density in this direction with a step function, centered 
around the peak of the halo. The width w r is deduced 
from a Gaussian fit to the radial profile of the halo den- 
sity obtained numerically. 

As shown in the Appendix [C] the assumptions intro- 
duced above lead to the approximate expressions 

V 2 = 1 + \ (h cl f?f?f? - h hh f* h fi h f* h ) ■ (26) 

Here we have introduced the function 



Y=1.02 



;l/bb 



(27) 



and the normalized bin widths are u^{ bb = L 7 _ f 1 2er 



cl/bb 



radial direc- 

l/bb 



in the z, t directions, while in the 
tion the limited density manifests itself by Ur 
min(L r , w r ) j 2a^^ hh . The limiting behavior of / is 

1 — it 2 /3 for small bins u much narrower than the cor- 
relations, and (l/«)yV/2 for large bins u — > oo. The 
above paramctrization is convenient, though elaborate. 
However, when bin widths tend to infinity, we obtain a 
particularly useful expression 




V [k J Q ] 



FIG. 3. The number squeezing parameter rj 1 in the dense 
7 = 1.02 case as a function of bin volume (log scale) calculated 
at three different scattering times (in units of t C oi). The solid 
lines are from a numerical solution of the full model (|5I4I) . The 
dashed lines are predictions of the Gaussian correlation model 
(|26[) . The dash-dotted horizontal line denotes the shot-noise 
level rf = 1. 



1.02 case 



rf = 1 + (27r)i^ [h^afafaf - W b W bb ] . (28) 

Hence, in general, in the multi-mode system, number 
squeezing depends on the number of particles in the bins 
n, but also depends on the correlation amplitudes hi and 
the widths of the correlation functions af^ hh . The situa- 
tion simplifies dramatically when the bins are very small. 
In this case we get 

V 2 = 1 + \ (ha - h hh ) , (29) 
and the quantum state is number-squeezed as long as 

h c l < hbb- 

To place this in context, an additional important re- 
mark is in order. If we consider a pure two-mode pair 
production model governed by the Hamiltonian H = 
ab + aW we obtain simply i] 2 = 0. Therefore we see 
that in a multi-mode system one cannot use the intu- 
itions from the two-mode model to predict such quanti- 
ties as the number-squeezing parameter, even when the 
bins being considered are very small. 



V. SIMULATION RESULTS AND ANALYSIS 

In this Section we perform a systematic study of the 
number squeezing parameter in condensate collisions and 
identify the physical phenomena, which have the largest 
impact on i] 2 . 



We begin with the case of a dense mother cloud, with 
7 = 1.02, so that the bosonic enhancement comes into 
play at some stage during the collision. In Fig. [5] we plot 
rj 2 as a function of time for three different bin volumes 
V (for details how the bins are chosen, refer to Appendix 
[A"]) . The solid curves, which are a result of the simula- 
tion of the full positive-P Equations (JSJ) are compared 
with the analytical prediction (|26|) based on correlation 
properties. The latter requires both the height and the 
widths of the second order correlation function as input, 
and these are extracted from the simulation. The Figure 
shows very good agreement between that model and a 
"direct" evaluation of the number squeezing by counting 
the number of particles in the two bins. Note that apart 
from early times, rj 2 > 1 so opposite bins do not reveal 
number squeezing. 

To check to what extent the number squeezing is a 
result of a particular choice of bins, in Fig. [3] we plot 
rj 2 as a function of the bin volume V at three times. 
Again, the agreement between the model (|2"o| and the 
simulation is very good. This Figure confirms that the 
absence of number squeezing in the scattering halo is a 
general feature regardless of binning details. 

In the past, it has been conjectured that the quantum 
correlations in the halo in such experiments are degraded 
by interactions with the mean field, or due to the time- 
variation fo the source cloud. To identify which process 
is responsible for such dramatic loss of number squeez- 
ing at later times, we first compare the results of the full 
positivc-P method with a maximally reduced Bogoliubov 
Method (RBM). The evolution of the colliding conden- 
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FIG. 4. The number squeezing parameter rf in the dense 7 = 
1.02 case as a function of time for three different bin volumes 
(in feo units), as predicted by the reduced Bogoliubov model 
(RBM) (|15I30|) . The dash-dotted horizontal line denotes the 
shot-noise level rj 2 — 1. 
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FIG. 5. The number squeezing parameter rf in the dense 7 = 
1.02 case as a function of time for three different bin volumes 
given in kg units. The solid lines are from a numerical solution 
of the full model H5I4|) . The dashed lines are predictions of 
the model without Bose enhancement (|8I9[) . The dash-dotted 
horizontal line denotes the shot-noise level n = 1. 



sates is simplified to a counter-propagation of the two 
initial clouds with fixed shape, and additionally, in equa- 
tions ([5]) we include only the pair production term, so it 
simplifies to (fl"5j) . Both the mean field self-interaction of 
the BECs and its impact on the scattered particles are 
neglected. Free kinetic dispersion is also suppressed, so 
that the equations of motion for the condensate field arc 

iWtM^t) = ^-(T2id x - fco)0±(x,t). (30) 

The condensates do not spread and the scattering pro- 
cess is maximally simplified. Figure U shows the number 
squeezing parameter as a function of time as predicted 
by the RBM. Although rj 2 does not grow as strongly as in 
Fig. [2j still the atom-difference fluctuations surpass the 
shot-noise limit. Therefore, it is neither the mean-field 
interaction nor the spreading of the BECs that have the 
major impact on the number squeezing parameter. 

Next, we simulate the condensate collision using the 
numerical method which a priori excludes the bosonic 
enhancement, introduced in Sec. (|IID|) . sec Eq. ©. In 
Fig. [5] we compare the results obtained in the full positive- 
P simulation and the non-enhanced method. We plot 
rj 2 as a function of time for three different bin volumes. 
Although the growth of the number squeezing parameter 
is less violent in the absence of bosonic enhancement, still 
rj 2 is above the shot-noise limit. As expected, for short 
times t < 0.1 the outcomes of the two methods agree very 
well - the system is still in the spontaneous regime. 

Finally, in Fig. [5] we draw the peak value of the back- 
to-back correlation function, namely the hj,b defined in 
Eq. l[25]). As indicated by equations {26} and {28j), the 
number squeezing is more likely to occur for high h b b, 
the widths of the correlation functions being the other 
factor. We see that in all three methods (full positive- 



P, RBM, and non-enhanced) give hu, < h c i = 1 at long 
times, a limit below which (|2T)|) indicates that small bins 
can never be number-squeezed. 
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FIG. 6. The peak height of the back-to-back correlation func- 
tion hbb as a function of time for the dense case 7 = 1.02. The 
solid line is from a numerical solution of the full model (|5I4|| . 
the dashed line from the model without Bose enhancement 
(|8l9p . while the dotted line comes from the reduced Bogoli- 
ubov model (RBM) (|15I30|I . The dot-dashed line shows the 
border value of /ibb = 1 when the back-to-back and collinear 
peaks are equal, and small-bin number-squeezing disappears. 
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FIG. 7. The number squeezing parameter rf in the dilute 
7 = 0.24 case as a function of time for three different bin vol- 
umes given in fcjj units. The solid lines come from a numerical 
solution of the full model (|5I4[) . The dashed lines are predic- 
tions of the Gaussian correlation model (|26p . The dash-dotted 
horizontal line denotes the shot-noise level rf = 1. 



B. 7 = 0.24 case 

Here we investigate the number squeezing parameter 
in the dilute case, when 7 = 0.24 and the bosonic en- 
hancement is absent. First, in Fig. [7] we plot if as a 
function of time resulting from the positive-P method 
§5§ and from the Gaussian correlation model (|2l)|) . The 
agreement is satisfactory, although the predictions of the 
model are very noisy. This is a result of the small number 
of scattered atoms. When the signal is low, it is difficult 
to extract the widths and peak values of the correlation 
functions that enter into the model. Nevertheless, we 
observe a major difference between the dense and the di- 
lute case. In the latter, the bins are number-squeezed, 
irrespectively of their volume and the time. 

We confirm that the system is indeed in the sponta- 
neous regime, by comparing in Fig.[5]the number squeez- 
ing parameter as predicted by equations ([5]) and the non- 
enhanced method ©. We do not observe any major dis- 
crepancy between these two results, and conclude that 
the system is indeed in the low-gain regime. 

To understand why the number squeezing is present 
in the 7 = 0.24 case, in Fig. [9] we plot the peak of the 
back-to-back correlation function as a function of time. 
We see, that it is far from reaching the border value of 
hbb = h c i = 1. Also, we compare this value with the 
one predicted by the non-enhanced method ((9]) and find 
excellent agreement. The number squeezing parameter 
depends not only on the peak values of the correlation 
functions, but also on their widths. In Fig.[TU]we compare 
the widths of the correlation function for 7 = 1.02 and 
7 = 0.24 as a function of time. We notice that in both 
situations the back-to-back and co- linear widths are com- 
parable to each other. In the dense case, the back-to-back 
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FIG. 8. The number squeezing parameter rj 1 in the dilute 
7 = 0.24 case as a function of time for three different bin 
volumes given in ko units. The solid lines result from a nu- 
merical solution of the full model (|5I4[) . The dashed lines 
are predictions of the model without Bose enhancement (18191) . 
The dash-dotted horizontal line denotes the shot-noise level 
„ a = l. 

h hh ^0.24 




FIG. 9. The peak height of the back-to-back correlation func- 
tion hbb as a function of time for the dilute case 7 = 0.24. 
The solid line is from a numerical solution of the full model 
(|5I4[) , the dashed line from the model without Bose enhance- 
ment (|8l9p . The dot-dashed line denotes the border value of 
hbb = 1 when the back-to-back and collinear peaks are equal. 



widths are slightly larger than the co-linear, which should 
favor rj 2 < 1 for large bins, as indicated by Eq. (|2"5)) . 

The lack of number squeezing for any bin size for long 
times of the 7 = 1.02 case, must be therefore attributed 
to the drop of the peak height hbb- This drop is due 
to the non-monochromatic nature of the parent BECs. 
Their momentum spread leads to a non-zero width of 
the back-to-back correlation function, which in turn re- 
sults in scattering into not exactly opposite momentum 
modes. The pair of atoms can end up in non-opposite 
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FIG. 10. The widths of the co-linear (left column) and back- 
to-back (right column) correlation functions for the dense (top 
row) and dilute case (bottom row) as a function of time. Fit- 
ted to the numerical solution of the full model (|5I4[1 . The solid 
black line is the width in the axial direction <r z , the dotted 
red line is in the tangential direction to the halo center at, 
while the dashed blue line is in the radial direction oy. The 
dot-dashed violet line is the Gaussian fitted half-width of the 
halo density w r , which narrows in time due to energy-time 
uncertainty principle. All widths are given in fco units. 



bins. Nevertheless, when the number of scattered atoms 
is low - as in the dilute gas case or at early times in the 
dense gas case - there is a high probability of finding a 
single or a few correlated pairs in the opposite regions, 
which is related to a high value of hbb- Crucially, the 
probability of having some uncorrelated pairs in the op- 
posite bin is low, because the neighboring bins are mostly 
empty. However, when the number of scattered particles 
grows, the chance that neighboring bins are empty be- 
comes small, and uncorrelated atoms spill over into the 
opposite bin. In this way, the rj 2 fluctuations grow, since 
there is a significant amount of uncorrelated pairs in the 
opposite bins. Figure [1] which is an outcome of a single 
collision in the dense gas case, is an illustration of this 
scenario. There are some clearly visible regions, where 
the atoms form a single large speckle, while on the other 
side of the halo two distinct speckles are present. 



i.e. (10/7) p. In Fig. [II] we plot rj 2 as a function of time 
for three different bin volumes. In all cases, the number 
squeezing is near perfect (if rts 0.03), despite a huge total 
number of scattered atoms ~ 10 7 . The residual slightly 
nonzero value of rj 2 is induced by the presence of a sea of 
short-lived particles from virtual scattering events. The 
RBM, which does not include this non-resonant effect, 



gives rj 
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FIG. 11. The number squeezing parameter n in the dense 
7 = 1.02 case as a function of time for several different bin 
volumes given in fcjj units. Here, the BECs are replaced with 
the plane-waves. The solid lines come from a numerical solu- 
tion of the full model (|5I4[) . while the dashed green line from 
a reduced Bogoliubov model (RBM) simulation (|15I30I) . The 
dash-dotted horizontal lines denote the shot-noise and zero 
levels. 

We also plot the peak height of the back-to-back func- 
tion as a function of time, see Fig. [T3] At long times it 
becomes indistinguishable from the border value of unity, 
but at this stage the number of atoms in the halo and the 
bin occupation n is very large. However, from the model 
([2TH one sees that the minimal, fully squeezed value of 
i] 2 = corresponds to hhb = 1 + 2/ri, which is extremely 
close to unity, so that this remains consistent with the 
observed strong number squeezing in Fig. 1111 



VI. CONCLUSIONS 



C. Collision of two plane- waves 

To confirm the conjecture formulated in the previ- 
ous paragraph, we simulate a collision of two monochro- 
matic plane waves. We use the same parameters as 
in the 7 = 1.02 dense case, but replace the initial 
BECs with monochromatic peaks in momentum space at 
k = (0, 0, ±fco) (and replace \(f>\ 2 in the evolution equation 
(HJ) with the mean density). This mean density is cho- 
sen equal to the peak density in the usual 7 = 1.02 case, 



We have performed a systematic study of the number 
squeezing parameter between two regions of opposite mo- 
menta in the scattering halo formed during collisions of 
two BECs. We have shown that the number squeezing 
depends strongly on the bin size, density of the mother 
clouds, mean-field interactions, and above all on the spec- 
tral purity of the mother clouds. In the dilute case, the 
number squeezing is evident, since the number of scat- 
tered pairs is low. Therefore, once a single atom is de- 
tected at momentum k, there is a high probability of 
finding one (and just one) at k' w — k. However, having 
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FIG. 12. The peak height of the back-to-back correlation 
function hbb as a function of time for the dense case 7 = 
1.02 and two colliding plane-waves. The solid line is from a 
numerical solution of the full model (|5I4|I . while dot-dashed 
line denotes the border value of hhb = 1 when the back-to- 
back and collinear peaks are equal. 

only a single - or a few - correlated pairs is unattrac- 
tive from a quantum intcrfcromctry point of view. The 
way to increase the number of scattered atoms is to use 
a dense condensate. In this case the system can enter a 
regime of bosonic amplification during the collision. 

However, as we demonstrated, there is a crucial draw- 
back of this approach. Namely, when the amount of 
scattered atoms is high, rj 2 grows rapidly. By perform- 
ing additional simulations with plane-waves instead of 
Gaussian condensates, we have related this behavior to 
the non-monochromaticity of the colliding BECs. As the 
sources are not point-like in momentum space and the 
back-to-back correlation function has a non-zero range, 
the scattered pairs are not perfectly aligned. This in turn 
results in imperfect correlation between opposite bins, 
which is additionally amplified by the increased num- 
ber of stray atoms that enter them from neighbouring 
regions. The non-monochromatic nature of the source 
clouds appears then to be the effect that underlies most of 
the degradation of pair correlations and number squeez- 
ing in the halo. Our results have potential application 
for the setup and analysis of existing and future experi- 
ments for the production of correlated atomic pairs from 
ultracold atom sources, as these are almost always multi- 
mode, i.e. non-monochromatic to an appreciable degree. 
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Appendix A: Division of the halo into bins 

To calculate the number-squeezing parameter we fol- 
low a similar procedure to that used in recent experi- 
ments [l|, We take an annular "washer-shaped" re- 
gion matched to the position and width of the halo, that 
excludes regions near the condensates, and divide it into 
zones. The matched dimensions of the washer-shaped 
region in the various calculations are shown in Table HI 

This annular region is then divided into equi-sized bins 
by making a series of equally-spaced cuts in the axial (z), 
radial, and tangential directions. Then number of cuts is 
varied to achieve a gradation from large bin sizes cover- 
ing an angle of it/2 in the k x — k y plane (8 bins in total, 4 
pairs) , to bin sizes comparable with the dimensions of a 
single correlation volume, then on to small bins the size 
of a single computational lattice point. The progression 
of bin dimensions is that we first reduce the radial bin 
size to approximately the correlation length ay , then the 
tangential and axial (z) sizes to their correlation lengths 
in that order. Then we repeat reductions in the same or- 
der down to single-computational-latticc-point volumes. 
One example is shown in Fig. [13] for the first 7 = 1.02 
dense case in the table. 

Number squeezing between opposite bins is calculated, 
then averaged over all pairs to obtain the displayed values 
of rj 2 . One subtlety should be mentioned: at small bin 
sizes, mean occupation can vary appreciably as a result 
of mismatches between bin shapes and the positions of 
discrete lattice points on the computational lattice; some 
small bins miss all lattice points altogether. Such an ef- 
fect does not appear experimentally, so to exclude poten- 
tially disruptive contributions from atypically discretized 
bins we exclude some bin pairs from consideration. The 
excluded bins are those for which either bin has an en- 
semble average occupation less than half or more than 
twice the mean bin occupation (when averaged over all 
bins). 
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FIG. 13. The bin widths (black solid line - axial, dotted red 
line - tangential, dashed blue line - radial) as a function of 
the total bin volume (dense case 7 = 1.02) used in the full 
model (|5I4[) simulations. 



Appendix B: Calcluation of halo correlations 



weights the contribution to the correlations proportion- 
ally to the product of the densities at the two points k 
and k+ Ak, which is approximately the local halo density 
squared. Thus, it takes into account primarily the most 
relevant, dense, part of the halo. Then, Gaussian fits are 
made to obtain the peak values h's and the widths cr's. 
During the fitting, points with excessive statistical uncer- 
tainty, or for distances Akj so large that the correlation 
function g^ begins to rise, are excluded. 
Appendix C: Gaussian model for halo correlations 

We use Gaussian g^> functions as in Equations (j2"4")l 
and (|25p and approximate the halo by a step function 

in the radial direction, k r = ^Jk^, + k^ . The density is 

Po when \k r — kh\ < Wh an d zero otherwise. Here, kh 
is the halo mean radius, and Wh the radial halo half- 
width taken to be \J n/2 times the standard deviation of 
a Gaussian fit to the true radial profile of the halo density. 
We also take the bins to be centered radially at kh so the 
effective integration range in the radial direction extends 
from —q r to +q r , where q r = min(L r /2, Wh). Therefore, 
the integrals of the correlation functions read 



The correlation properties are averaged in a similar 
way to that performed in experiments [l], 0]. That is, 
we calculate the mean correlation function in a region 
IZj as a function of inter-particle distance Afc in various 
directions j in the following way: 



L z /2 L t /2 



G aa /ab = Po J dk z dk' z J dk t dk' t J dk r dk' r g^j hh (k,k'), 

-L z /2 -L t /2 -q r 

which - using Equations (|T8|) and (fT9| - gives 



_ (2) / TC . d 3 k(^t ( k ) ( k + Ak)S(k)6(k + Afc)) 

■^ cc / bb( '~ J n ,d 3 kn(k)n(k + Ak) 

(Bl) 

This region is similar to the annular region used to calcu- 
late number-squeezing (as shown in TablcUand explained 
above in Appendix [XJ for the j = z direction, but a re- 
duced volume for the other directions j = t, r, where the 
reduction consists of additional restriction to within fco/3 
of the x and y axes. However, the washer shape is much 
wider radially, i.e. it has the same average radius as in 
the table, but a radial width of 0.9fco- 

The averaging method used here is as employed in ex- 
periments because it is convenient for analysis of data 
consisting of detected particle positions. It effectively 



l /ab 



i + ^/bb/K 1/bb )/K 1/bb )/«. 1/bb ) 



Here 1 



/(«) = i 



—u erf (it\/2) — — ^1 — e~ 



is a function of the normalized bin widths 

cl/bb _ L x ,t ..cl/bb 1r 



L z,t 



cl/bb 



U 



cl/bb ' 



This expression is inserted into Eq. (|22j) to obtain the 
estimates of rf 1 on the basis of correlation and density 
measurements. 
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